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We study the nonequilibrium aging dynamics in a system of quasi-hard spheres at large density 
by means of computer simulations. We find that, after a sudden quench to large density, the relax- 
ation time initially increases exponentially with the age of the system. After a surprisingly large 
crossover time, the system enters the asymptotic aging regime characterized by a linear increase of 
the relaxation time with age. In this aging regime, single particle motion is strongly non-Fickian, 
with a mean-squared displacement increasing subdiffusively, associated to broad, non-Gaussian tails 
in the distribution of particle displacements. We find that the system ages through temporally inter- 
mittent relaxation events, and a detailed finite size analysis of these collective dynamic fluctuations 
reveals that these events are not spanning the entire system, but remain spatially localized. 
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I. INTRODUCTION 

Aging refers to the slow evolution with time of the 
physical properties of a disordered material suddenly 
quenched into a glass phase This might refer to the 
evolution of the density of a polymer glass, of the di- 
electric response in an organic liquid, or the height of 
a gently shaken sandpile. Aging is more easily detected 
by focusing on the dynamics of the system, for instance 
by measuring correlation or response functions. One can 
think of measuring the decay of density fluctuations us- 
ing light scattering techniques in colloidal glasses, or the 
timc-depcndcnt magnetic response in a spin glass. After 
several decades of aging studies in glassy materials, these 
generic features are very well documented [H-Q . 

Much less is known, however, about the microscopic 
mechanisms involved during aging, and how these evolve 
with time, for two main reasons. First, the important 
theoretical developments from the last decade about ag- 
ing dynamics mainly stemmed from 'mean- field' types of 
approaches, which provided detailed predictions about 
the behavior of averaged dynamic quantities, but very lit- 
tle information about microscopic motion Second, it 
is only relatively recently that the microscopic dynamics, 
of, say, molecules in a supercooled liquid, has been char- 
acterized in detail in equilibrium conditions above the 
glass transition and, by comparison, much less work 
has been done about the corresponding aging regime at 
lower temperatures. 

In this paper, we use numerical simulations to study 
the aging dynamics of a model system designed to un- 
derstand the behavior of dense suspensions of colloidal 
hard spheres. Hard spheres represent one of the simplest, 
and thus most studied, models to study the glass tran- 
sition. It is easily studied in simulations, and the model 
finds its experimental realization using well-controlled 
colloidal suspensions [11,0]. As opposed to molecular liq- 



uids, colloidal hard spheres can be studied at the particle 
scale using microscopy techniques 0, Q j while dynamic 
light scattering provides a convenient way to characterize 
their dynamics in great detail [l^ . It should be noted 
that while a relatively small set of experiments have been 
reported on the aging of colloidal hard spheres jol-fisj. a 
much broader literature exists on the out-of-equilibrium 
dynamics of a variety of more complex colloidal systems. 
Throughout this paper, we will in particular make refer- 
ence to findings for colloidal attractive gels [l3 - [T6j and to 
systems with soft repulsive interactions, such as closely 
packed soft spheres [17Hl9l | or Laponite clay platelets in- 
teracting via Coulomb repulsion |2C)| - [24| . 

In this work, we will be concerned with three main 
questions. 

(1) How does the structural relaxation slow down with 
the aging time? While numerical simulations of model 
glasses gcncrically find that the structural relaxation time 
increases roughly linearly (or sub-lincarly) with the sam- 
ple age psj . some light scattering experiments on gels 
and Laponite report an unexpected exponential growth 
of the relaxation time with sample age [3 [l^, 21 1, at 
least at short times. 

(2) How do particles move during the aging regime? 
Confocal microscopy experiments report that colloidal 
particles move very little in concentrated hard spheres, 
to the point that distinguishing thermal vibrations 
from gen uine relaxation becomes an experimental chal- 
lenge jl2l . [isj . In simulations as well, particles move very 
little, leading to a very slow, troically algebraic, decay 
of time correlation functions [25l - [27| . This is in stark 
contrast with several experimental reports of sharply de- 
caying time correlation functions in aging molecular liq- 
uids seeded with colloidal particles or in colloidal gels 
and Laponite, where compressed exponential decay and 
ballistic particle motion over large distances were re- 
ported [IHIli,!!!,!!!. 
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(3) How heterogeneous is the dynamics? Due to the 
key role played by dynamic heterogeneity in equilib- 
rium studies of the glass transition Q , similar signatures 
have been sought in aging materials. Optical and con- 
focal microscopy studies revealed the existence of rare 
small-scale relaxation events involving extended clusters 
of particles 0], whose size is modest and grows very lit- 
tle , or not at all [13, Ell , during aging. This agrees 
with numerical simulations on model glasses, where four- 
point dynamic susceptibility appear to grow at a very 
slow rate with the sample age [13, HH, suggesting that 
if collective displacements occur, they are correlated 
over relatively short length scales. This set of results 
is however in contrast with another series of observa- 
tions. Sudden relaxation events termed (rather dramati- 
cally) 'earthquakes' [11] or 'avalanches' |32l . |33| were re- 
ported in numerical simulations of Lcnnard- Jones glasses 
and these were suggested to be spanning the entire sys- 
tem. Similarly, highly intermittent dynamic fluctuations 
were experimentally recorded in several aging systems 
from polymer glasses [s^ to colloidal gels [ij. El, [13 
and soft spheres [ll], the latter typically involving bal- 
listic particle motions correlated over very large dis- 
tances [H El El El- 

The paper is organized as follows. In Sec.|TTl we present 
the numerical model and techniques used in the present 
study. In Sec. IIIIl we describe the average aging dy- 
namics for the evolution of the energy and relaxation 
timescalc. In Sec. IIVI wc focus on the subdiffusive and 
heterogeneous dynamics at the single particle level. In 
Sec.|V]we present results concerning the intermittent col- 
lective dynamics of the system. In Sec. lVII we discuss our 
results and conclude the paper with some perspective for 
future research. 



II. NUMERICAL MODEL AND TECHNIQUES 

We perform molecular dynamics simulation of dense 
systems of strongly repulsive particles interacting with a 
very steep pair potential designed to model the behavior 
of colloidal hard spheres [11] : 

where e is an energy scale, rij — jr^ — rj|, being the 
position of particle i, and atj = (ct; -f- crj)/2, where ai 
represents the diameter of particle i. To prevent crystal- 
lization occurring in dense systems of hard spheres, we 
introduce a size polydispersity and draw the particle di- 
ameters from a flat size distribution, ai/a £ [0.8, 1.2], so 
that the average diameter is o7 = cr, and the polydisper- 
sity S is given by 



6 = -1^11.5%. (2) 
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To study the dynamics of the system we solve Newton's 
equations for a system composed of N particles, using a 
velocity Verlet algorithm |36[ in a cubic box of linear size 
L. We use periodic boundary conditions in the three 
directions of space. We have studied two system sizes, 
N — 4000 and N = 500. When dealing with averages, 
wc shall report results for the largest system size, while 
a comparison between the two systems will allow us to 
perform a detailed comparison of the dynamic fluctua- 
tions occurring during the aging dynamics as a function 
of system size. We perform simulations in the NVT en- 
semble, and we control the temperature by rescaling ve- 
locities every 100 molecular dynamics timesteps in order 
to maintain the kinetic energy to the desired constant 
value. 

For the inverse power law potential in Eq. ([T]), density 
and temperature cannot be controlled independently, as 
rescaling the density by a factor A is equivalent to rescal- 
ing the temperature by a factor A~^/^^. Thus we fix the 
temperature and energy scales, ksT = e = 1/3, and sim- 
ply vary the density, p = N/L^, of the system [11]. To 
ease the comparison with hard sphere experiments, we 
express our results using the volume fraction, ip, rather 
than density p , where 

i 

In the following, we express length scales in units of a, 
and timescales in units of tq = a^mlt. where m is 
the mass of the particles. We use a time discretization 
Ai = O.OItq, which ensures a proper integration of the 
equations of motion. 

We shall study the aging dynamics of samples in the 
range = 0.553 — 0.662. To obtain reproducible results, 
we need to produce fully disordered initial states. To 
this end, we first equilibrate the system at very low vol- 
ume fraction, (p ~ 0.14. Wc then compress the system 
very rapidly to the desired final volume fraction. The 
compression is done in small successive steps in order to 
avoid large overlaps leading to very large repulsive forces. 
The age t„ of the system is counted from the time when 
the system reaches the final volume fraction. To increase 
the statistical significance of our results, we have per- 
formed 5 independent runs at each volume fraction with 
N = 4000 particles, starting from independent configura- 
tions compressed from the fluid at = 0.14. During the 
course of the simulations we found no sign of incipient 
crystallization in the system, while crystallization was a 
major obstacle in an initial set of studies using a smaller 
polydispersity, (5 w 6 %, for which we do not present 
results. 

We have previously studied the equilibrium dynamics 
of this system [ll] . We found that the dynamics slows 
down considerably when ip increases above 1^ « 0.50. 
Fitting the increase of the relaxation to a power-law di- 
vergence at some critical volume fraction (/3c, we obtained 
if>c ~ 0.592. This fit then locates the (ap parent) mode- 
coupling singularity for this system |37{ . which should 
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FIG. 1: Aging of the potential energy density, Eq. ((4]), for 
= 0.647 and (p — 0.662. Data are equally well fitted with 
algebraic, Ai x t"" + _Bi, or logarithmic, A2 x (Intw)*" + -B2, 
time dependences for waiting times larger than a surprisingly 
large crossover timescale, Ta, indicated by arrows. 

serve as a reference volume fraction for the aging studies 
below, since it becomes very difficult to reach, thermal 
equilibrium within the accessible numerical timescales 
when if increases beyond ipc- 



III. TOWARDS THE ASYMPTOTIC AGING 
REGIME 

We now present the results of our study, starting from 
the time evolution of the simplest quantity we can mon- 
itor in our study, the potential energy density, defined 
as 

Miw) = (^EE^(^».))' (4) 

\ i=i j>i I 

where the brackets stand for an average over independent 
initial configurations. Wc present representative results 
for the time evolution of epit^) at large densities in Fig.[T] 
The central observation from this figure is that, at 
these large densities, the potential energy keeps evolving 
over the 6 decades of the simulation and never reaches an 
asymptotic plateau. This simply confirms that thermal 
equilibrium cannot be reached at these volume fractions, 
because the equilibrium relaxation time is much larger 
than the maximum time accessible to our simulations. 
More in detail, we also observe that the energy evolves 
initially quite rapidly, decreasing by a factor of about 2 
when time increases from ^ 1 to = 10^, while it 
evolves only by a few percent when time further increases 



by two additional decades. Thus, we clearly see the effect 
of 'aging', since dynamical evolution slows down consid- 
erably as the age of the system increases, as is well-known 
from decades of experimental aging studies in many dif- 
ferent materials. 

To describe the slow evolution of the potential en- 
ergy, we fitted our data to both a power law decay, 
ep(tw) ~ ^i^w" + -^ii ^^"i ^ logarithmically slow evo- 
lution, ep{t^) = (Iniw)'' + S2. Wc find that both 
fits give an equally good description of the final, slow 
evolution of the potential energy. The fit parameters, 
a ?» 0.2 — 0.3 and b « 0.1 confirm that the evolution of 
the energy is indeed extremely slow. 

We also find that the fits only hold when is very 
large: the first few decades of the simulation, correspond- 
ing to a faster evolution of ep(t„) are not well described 
by this asymptotic slow decay. This implies that a rel- 
atively long time is needed for the system to enter the 
asymptotic aging regime. Moreover, we find that the 
time Ta required to enter the asymptotic regime increases 
when density increases and the system is quenched deeper 
into the glassy state. This implies that considerable care 
must be taken when performing data analysis of the aging 
regime, since it already takes a large part of the simula- 
tion simply to enter the final regime, indicated by the 
arrows in Fig. [TJ It is likely that the first regime corre- 
sponds to faster processes where the largest overlaps cre- 
ated during the compression are removed, producing heat 
which is then removed by our thermostatting procedure. 
Thus, no 'universal' characteristics are to be expected in 
this time regime, which might well depend quite strongly 
on the details of the preparation procedure, or the cho- 
sen thermostat. Note that Ta ~ 10^ ^ tq for the data 
presented in Fig. [1] meaning that scaling behavior might 
only be observed in the demanding limit of 

tw > T„ > To. (5) 

It is crucial to recognize the existence of two distinct 
aging regimes in order to analyze properly the scaling 
properties of dynamic functions in the asymptotic ag- 
ing regime, as we find that these two distinct regimes 
in fact affect most of the measurements we made in our 
simulations. In Fig. [5] we show the behavior of the self- 
intermediate scattering function following a quench at 
if = 0.647. This two-time quantity is defined as 

/.(iw,r) = ^lge^'l-[r.(*w+r)-r,(*,)]^ (g) 

and we perform measurements at q ~ 7.8, close to the 
first peak in the static structure factor. As is well-known, 
two-time quantities reveal the non-stationary evolution of 
the system in a very direct manner. In agreement with 
previous work (25l . \2l\ , we find that the decay with the 
delay time r of /^(fwj t) occurs in a two-step manner, the 
slow decay being a strong growing function of the waiting 
time tw. The faster initial decay is much less sensitive to 
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FIG. 2: Self- intermediate scattering function, Eq. ([6]), for 
= 0.647 in a log-log representation. The long-time decay 
is fitted with two power law decays for tw < Ta, a single one 
for > Ta, as shown with full lines. The exponent for the 
latter is nearly constant, ~ —0.15. The horizontal dashed line 
indicates fs ~ 1/e. 



FIG. 3: Evolution of the relaxation time, Ta{tm), for different 
volume fractions. The system at = 0.553 reaches equi- 
librium (horizontal line). A simple aging regime, Tc ~ tw, 
indicated by the full lines, is reached for larger densities, af- 
ter a transient regime where Ta grows exponentially with t^, 
indicated with a dashed line which ends for ~ Ta- 



time, 



the waiting time, and corresponds physically to particle 
vibrations in a nearly frozen amorphous structure. 

The existence of two distinct aging regimes is obvious 
from these data, since the slow decay of fs clearly shows a 
crossover for time delay r such that r-t-t^ ~ ^a, as shown 
in Fig. [21 while the data become simpler to describe when 

> Ta as the entire decay then takes place in the asymp- 
totic long-time regime defined by Eq. ([S])- As found in 
other systems [l^ , we find that the long-time decay of the 
self-intermediate scattering function is well described, in 
the aging regime, by a power law, fsit^^r) ^ t~", as 
indicated in Fig. [5] by the full lines. As noted before, 
such an algebraic decay is in contrast with equilibrium 
measurements at lower volume fraction which typically 
display stretched exponential relaxations. Stretched ex- 
ponential decays are only seen for shallow quenches when 
a crossover towards thermal equilibrium is possible. The 
small value of the exponent, a ~ 0.15 means that relax- 
ation occurs over a very broad time window, and that it 
is not possible to define a 'typical' relaxation time for the 
system. 

Nevertheless, to quantify the aging of the dynamics, 
we follow the common practice and define an a-rclaxation 
time Ta{t^) from the time decay of fg as /s(iw, t^) = 1/e, 
as indicated by the horizontal dashed line in Fig. [2] In 
Fig. [31 we report the evolution of Ta(tw) with waiting 
time for all volume fractions investigated in this work, 
from (fi = 0.553 to (p = 0.662. Again, we observe two 
distinct regimes for the evolution of Ta{tw). While for 

< Ta-, Ta sccms to incrcasc exponentially with waiting 



exp(ctw), <C Ta, 



(7) 



with c a numerical constant, its growth is better de- 
scribed by an algebraic law at long times, > Ta, which 
is well compatible with a so-called 'simple aging' behav- 
ior 0, 



> Ta- 



(8) 



This simple aging behavior is shown with full lines in 
Fig. [3l In agreement with the observations in Fig. [1] we 
find that it takes an increasingly long time for the system 
to enter the asymptotic (simple) aging regime when ip in- 
creases. For the larger (p studied, (p = 0.662, the simple 
aging regime is not reached during the course of our sim- 
ulations, and the system appears to be in the crossover 
between ([7]) and ([8]), and the system seems to undergo an 
effective 'super-aging', i.e. its relaxation time increases 
faster than its age. Finally, the opposite 'sub-aging' be- 
havior is found when volume fraction is not too large, for 
instance ip = 0.592, because the system is crossing over 
towards thermal equilibrium, such that should satu- 
rate at long waiting times to the equilibrium relaxation 
time. These observations show that a complex behavior 
of Ta, as often reported in numerical works 25, 38| and 
experiments on colloidal gels and Laponite IJ, [20, 



may be due to the occurrence of multiple crossovers which 
are highly sensitive to volume fraction. 

These observations are also in agreement with the com- 
mon observation of sub-aging behavior in aging molecu- 
lar liquids, for which experiments are traditionally per- 
formed not very far below the glass temperature [l|, [40| . 
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FIG. 4: Mean-squared displacements, Eq. at (f — 0.647 
with symbols as in Fig. [21 The dashed line indicates difTusive 
behavior. Inset: log-log plot of the data for > Ta and large 
r with subdiffusive fits, Eq. (jlOp . indicated with full lines. 



Exponential growth of the relaxation time sometimes 
(but not always) followed by algebraic growth, has been 
reported in a few experiments on colloidal glasses or gels 
as well [H, [13, Our results thus suggest that this 
exponential growth might well be a transient behavior, 
which can persist, however, over a very large time win- 
dow, in particular for very deep quenches. Our results 
also highlight the difficulty in analyzing the scaling prop- 
erties of two-time dynamic quantities in numerical sim- 
ulations, since the asymptotic aging regime is only ac- 
cessed after a time which can be very large, thus de- 
creasing the window where 'universal' aging properties 
can be probed. This might well explain some conflicting 
results and analysis reported in the literature regarding, 
e.g., the proper scaling behavior of the self intermediate 
scattering function in Lennard- Jones glasses [4l|, HH . 



IV. SUBDIFFUSIVE AND HETEROGENEOUS 
DYNAMICS 

The above findings that the energy decays very slowly, 
while time correlation functions decay in an algebraic 
manner indicate that there is actually very little dynam- 
ics taking place in the system. To confirm this idea, we 
measured the averaged mean-squared displacement. 



(9) 



In Fig. m we present the time evolution of the mean- 
squared displacements at volume fraction ip — 0.647 mea- 
sured at different waiting times, using both a standard 



log-log representation (inset) and a lin-log representation 
(main plot). As found in Fig. [2] for the self- intermediate 
scattering function, the dynamics proceeds again in a 
two-step fashion with a rapid, ballistic regime at short- 
time, corresponding to vibrations of the particles in a 
frozen amorphous structure, followed by a much slower, 
waiting time dependent structural relaxation, which be- 
comes slower when t„ increases. This implies that the ag- 
ing of the system corresponds, at the microscopic scale, 
to a dramatic slowing down of single particle displace- 
ments, as found in experiments [TTI.TT2l.lT7j . 

It is significant that in order to represent graphically 
the behavior of the mean-squared displacements, we had 
to use a linear scale for its amplitude, with a range which 
remains smaller than Ar^ =0.5. This implies that over 
the entire duration of the simulation, particles, on av- 
erage, move a distance which is smaller than the mean 
particle diameter: despite the non-trivial aging dynam- 
ics we discuss, it should be clear that the structure of 
the material remains in fact almost frozen as soon as the 
system is quenched in the glassy phase, with essentially 
no large-scale dynamics taking place. 

The linear representation in Fig. 2] also confirms the ex- 
istence of the two distinct aging regimes discussed above, 
thus directly revealing the influence on the dynamics at 
the particle scale of the timescale Ta- Again, the data ob- 
tained for < Ta present a crossover for a delay r given 
by T -I- ^ Ta, as seen in Fig. |4l Thus, in the following 
we again focus on waiting times that are larger than r^, 
which, we believe, reflect better the 'universal', quench- 
independent nature of the aging regime in concentrated 
hard spheres. 

In the long time regime, the dashed line in Fig. 21 which 



corresponds to a diffusive growth, Ar 



is obviously 



an incorrect description of our data since the displace- 
ments seem to increase much more slowly than diffusively. 
As shown in the inset, a subdiffusive law describes our 
results very satisfactorily. 



Ar^{t^,T) r^T"", t^:>Ta, 



(10) 



where < 1 is an exponent characterizing the subdiffu- 
sion. 

We have measured v over a broad range of waiting 
times and volume fractions, and we present its evolution 
in Fig. [5] Apart from the low volume fractions where 
(diffusive) thermal equilibrium is reached rapidly, we sys- 
tematically find that the system obeys subdiffusive rather 
than diffusive behavior. At a given volume fraction, the 
system gets closer to equilibrium for larger i^, and cor- 
respondingly we find that the exponent ly increases with 
the age of the system, although it always remains very far 
from its equilibrium, diffusive value v = I. Note that an 
increasing v does not imply that particles move faster at 
large waiting times, as the prefactor in the power law ((TU]) 
is itself a decreasing function of t^, which implies that 
the total amplitude of the particle displacements indeed 
decreases with the age of the system. 

As volume fraction increases, the system is quenched 
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FIG. 5: Evolution of the subdiffusive exponent v for differ- 
ent and ip. The exponent is further from its equihbrium 
diffusive hmit u = 1 for smaller tw and larger ip. 



deeper and deeper into its glass phase, and deviations 
from diffusive behavior are correspondingly larger, which 
translates into a subdiffusive exponent v which gets 
smaller when ip increases, see Fig. [5] This simply means 
that particles move less and less when ip becomes larger, 
which is physically expected. 

While the mean-squared displacement quantifies the 
average dynamical behavior of the system, it does not 
convey much information on dynamic fiuctuations, which 
are recognized as an important feature of the single par- 
ticle dynamics in glassy materials. To quantify this 
'dynamical heterogeneity' during aging, we measure the 
probability distribution function of the single particle dis- 
placement, also known as the self-part of the van- Hove 
correlation function: 



1 



N 



\ 1=1 I 

(11) 

where Xi{t) represents the projection of Ti{t) along the x- 
axis. To improve the statistics, we use the isotropy of the 
system and further average Gs(a^, tw, t) along the three 
directions of space. For simplicity, we keep the notation 
Gs{x,t^,T) for the van- Hove function averaged over the 
x, y , and z directions. 

To gain deeper insight into the subdiffusive behavior 
described above, we present in Fig. [5] the typical shape 
and evolution of the van-Hove function measured for 
tf = 0.647, for a large waiting time, = 2451 in the 
asymptotic aging regime, and a delay r within the long- 
time subdiffusive regime. The data are presented in a 
lin-log scale, where Fickian behavior leading to a Gaus- 
sian shape of the particle distribution would appear as an 
inverted parabola, as shown by the dashed lines. These 
data arc highly reminiscent of the equilibrium findings 
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FIG. 6: Distribution of single particle displacements, Eq. 
for if = 0.647, tw ~ 2451, and various r indicated in the 
figure. While a Gaussian (dashed lines) fits the center of 
the distribution, the tails are broader and well described by 
Eq. (|12p . as shown by the full lines. 



that van-Hove functions are well-described at short dis- 
placements by a Gaussian form, with tails that are much 
broader than the Gaussian prediction 0, |1| . This implies 
that a small fraction of the particles move significantly 
farther than what is expected from a purely diffusive and 
homogeneous process. This is the most basic observation 
that the system is dynamically heterogeneous and can be 
described, as being composed of distinct families of 'slow' 
and 'fast' particles, a well-established distinctive feature 
of many disordered, glassy materials [3] . 

How should we describe the functional form of the 
tails of these distributions? At thermal equilibrium, it 
is now understood that the tails are well-described by 
an exponential (rather than Gaussian) decay, Gs{x,t) ^ 
exp(— |a;|/A), for r corresponding to the alpha-relaxation 
timescale [43] . The physical origin of the exponential, 
detailed in Refs. |43l - |45| . is due to the stochastic nature 
of the diffusion process in disordered materials, which is 
well described by the formalism of continuous time ran- 
dom walks (CTRW) [i^. In this description, particles in 
a dense supercooled liquid undergo a succession of long 
periods of vibrations within the cages formed by their 
neighbors, followed by rapid 'jumps'. Both the size of 
the jumps and, more importantly, the timescale separat- 
ing them arc statistically distributed quantities. The ex- 
istence of these statistical distributions directly accounts 
for the exponential form of the tails in the van- Hove func- 
tions 



43|. 



The natural extension of these considerations to the 
non-equilibrium aging regime studied in the present work 
is the a ging continuous time random walk (ACTRW) for- 
malism [431 J whose main features are those of the CTRW 
recalled above. The only difference lies in the functional 
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FIG. 7: Evolution of the exponent /3 obtained by fitting the 
tails of the distribution of particle displacements to Eq. (|12|l , 
for tw = 2451 and several r and ^p. 



asymptotic subdiffusive aging regime. We remark that 
these data have considerably more scatter than the data 
for the subdiffusive exponent i/, confirming the difficulty 
to obtain a reliable determination of /? from our data. 
Regardless of the noise, we note that /3 values are in the 
interval f3 G [1,2], as expected from Eq. (fT3| . suggest- 
ing that our description is physically sound. Moreover, 
the data at large density present a systematic evolution 
of /?, which increases towards 2 with t^. Assuming that 
Eq. (fT3| is correct, this would imply that i/ increases to- 
wards unity with t^, as indeed is observed in Fig. (5] The 
data at moderate volume fractions, ip < 0.60, are more 
difficult to interpret as (3 seems to decrease with t^, in 
this regime, while ly was found to increase. This could be 
due to the fact that for (p = 0.592 w ipc, the subdiffiisive 
aging regime is crossing over towards equilibrium during 
the simulation, and so the results could be a subtle com- 
bination of the exponential tail reported at intermediate 
times for equilibrium dynamics j43j , and the subdiffusive 
regime found deeper in the glass and described in the 
present study. 



form used for the distribution of times separating the 
jumps, which acquires 'fat', non-normalizable tails in the 
aging regime, in direct analogy with the trap model in- 
troduced by Bouchaud to describe aging dynamics in 
glasses [i^ . In this approach, the tails of the jump time 
distribution decay as ~ t~^~'^, where v < 1 is the 
subdifFusion exponent introduced in Eq. (|10|). while the 
tails of the self-part of the van-Hove function are not 
exponential anymore, but are instead asymptotically de- 
scribed by [4^ 



Gs{x,t„,T) - |x|-''exp(-(|x|/A)^) 



(12) 



where the exponent /? is also connected to v via the re- 
lation 



/3 = 2/(2-z.), 



(13) 



so that /? = 2 (Gaussian) is recovered when 1/ = 1 (Fick- 
ian diffusion). For subdiffusive processes, < < 1, one 
expects instead 1 < /? < 2. 

To test in detail the ACTRW picture, one should mea- 
sure the distribution of jump times V'(0, ^-^^d use it to di- 
rectly predict the mean-squared displacements and van- 
Hove functions [50{ . This is an interesting project, but is 
beyond the scope of the present study. Instead, we more 
simply use Eq. ((T^ as a theoretically motivated fitting 
formula for the tails of the Gs distributions. As shown 
in Fig. [51 such a fit describes the tails very well. We find 
a similarly good description for a broad range of volume 
fractions and timcscalcs. Unfortunately our data is not 
accurate enough to allow for a very precise determina- 
tion of all the fitting parameters. Thus, the following 
results should be discussed at a qualitative, rather than 
quantitative, level. 

In Fig. [7] we present the evolution of the exponent 
/? describing the tails of the van-Hove functions in the 



V. INTERMITTENT DYNAMIC 
FLUCTUATIONS 

Our study of single particle dynamics suggests that 
particles undergo very little dynamics, leading to a sub- 
diffusive growth of mean-squared displacements, associ- 
ated to a broad distribution of particle displacements. 
The picture of ACTRW which we used to describe our 
data is very similar in spirit to Bouchaud's trap model, 
and both suggest that the aging dynamics of concen- 
trated hard spheres is a temporally intermittent pro- 
cess [511 . It is the aim of this section to establish whether 



this picture is indeed correct. 

To this end, we must turn to collective dynamic ob- 
servables, and resolve the aging dynamics in space and 
time, and focus on /g, the instantaneous (unaveraged) 
value of fg. We present in Fig. [Hjja) the results of inde- 
pendent realizations of the dynamics at a given waiting 
time, = 2451, and volume fraction, ip — 0.637. For 
this system, containing N = 4000 particles, we observe 
small run-to-run fluctuations, which show that resolving 
the temporal evolution of the dynamics in a system of 
linear size L ~ 15(t is not sufficient to reveal significant 
dynamic fluctuations. 

Thus, we improve our spatial resolution and repeat this 
analysis for a smaller system with L 7.5a which con- 
tains N = 500 particles, see Fig.[5i;b). The data in Fig.[5| 
show that run-to-run fluctuations of fs become larger 
when N is smaller, as expected. However, we empha- 
size that the nature of the dynamic fluctuations seems to 
change qualitatively when the size is reduced. For a sin- 
gle realization of a quench with a small enough number 
of particles, the decay of the self-intermediate scattering 
function in the long-time regime is highly intermittent, 
as shown by sudden 'drops' of fs [26j. From one run to 
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FIG. 8: Instantaneous (unaveraged) intermediate scattering 
function, /s(tw = 2451, r) for the system at 95 = 0.637 and 
5 independent quenches, (a) measured in systems with A'^ = 
4000; (b) N = 500. While small run-to-run fluctuations are 
observed in (a), a wide variety of sudden decorrelation events 
are observed in (b). 



another, both the temporal location and the amplitude of 
these drops vary considerably, indicating that in each run 
the dynamics proceeds via temporally intermittent relax- 
ation events that mobilize a finite fraction of the whole 
system. However, the comparison with the fluctuations 
obtained with = 4000 suggests that the spatial exten- 
sion of these intermittent events docs not increase with 
the system size, but remains instead spatially localized. 

Note that this finding could be due to a spurious fi- 
nite size effect. An alternative interpretation of the data 
in Fig. [8] could be that sudden decorrelations in the sys- 
tem only occur when the system is too small, and these 
events are not observed in a bigger system because they 
disappear altogether. To disprove this interpretation, 
we reanalyze a single quench with N = 4000 particles, 
and further decompose the computation of the correla- 




FIG. 9: Top: intermediate scattering functions /s(tw = 
2451, r) measured in sub-systems containing A*' = 500 par- 
ticles within a larger = 4000 system. The variability be- 
tween sub-systems in a single realization demonstrates that 
intermittent decorrelation events are spatially localized. Bot- 
tom: snapshot of the A'^ = 4000 system showing the 5% most 
mobile particles over the time interval shown by the two ar- 
rows in the top panel (r — 1174 and 2789, respectively) are 
drawn with a larger size. They occupy essentially a small 
corner of the simulation box. 



tion function in 8 distinct sub-systems, each comprising 
500 particles and corresponding to a cubic sub-box of 
linear size L/2. In Fig. |9l we present the result of this 
analysis. Remarkably, we find that in a single quench, a 
small part of the system might indeed undergo the type 
of sudden rearrangement seen m N ~ 500 particles runs, 
but decorrelation events are not necessarily seen in other 
subsystems. This is visible in the top graph, where repre- 
sentative fs for various sub-boxes are shown, and in the 
bottom snapshot, where the 5% most mobile particles 
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FIG. 10: The four-point susceptibility X4, Eq- (El, for t„ = 
2541, — 0.637, and two system sizes, shows no system size 
dependence, and its modest peak value, X4 ~ 10 confirms the 
localized nature of decorrelation events. 



between the two times shown by the arrows in the main 
plot are depicted as large spheres. Clearly, individual re- 
arrangement events are confined to a fraction of the total 
volume (in that case a corner of the box), thus affecting 
significantly only the correlation function of the corre- 
sponding sub-box. We conclude that in a large system 
the aging dynamics occurs through an intermittent suc- 
cession of spatially localized decorrelation events, with 
no obvious 'confinement' effect imposed by the use of too 
small a system size, at least for the two system sizes used 
in this study. 

To confirm quantitatively both the absence of seri- 
ous finite size effects, and the spatially localized nature 
of decorrelation events, we have measured the variance 
of the spontaneous fluctuations of the self-intermediate 
scattering function, also known as a 'four-point dynamic 
susceptibility (30ll3]|: 



X4(iw,T) =N({fs (t„,T)) - (/,(t„,T))- 



(14) 



where /, represents the instantaneous value of fs- With 
this normalization, X4 should become independent of the 
system size in the thermodynamic limit, with an am- 
plitude that gives direct access to the typical number 
of particles relaxing in a correlated manner during ag- 
ing [H,!!!]. The data for N = 500 and TV = 4000 shown 
in Fig. [TU] are essentially the same, within our statistical 
accuracy, for both system sizes, and Xi reaches a value 
near xa ~ 10 for timescales corresponding to the slow re- 
laxation of the correlation function, confirming the spa- 
tially localized nature of relaxation events. 



VI. DISCUSSION 

We have studied numerically the aging dynamics of a 
concentrated system of nearly hard spheres over a broad 
range of volume fractions, and a large time window. We 
have addressed the three central questions presented in 
the introduction, concerning the evolution of the dynam- 
ics, single particle motion, and the heterogeneous nature 
of the dynamics. 

We found that over the 6 decades of aging times cov- 
ered by our simulations, the structural relaxation con- 
tinuously slows down, as expected for any aging system, 
but we showed that an asymptotic aging regime is only 
accessed after a timescale r^, which can become quite 
large for dense systems. Interestingly, we found that the 
relaxation time Ta{t^) first increases exponentially, as ob- 
served in some experiments on colloidal gels and repulsive 
Laponite platelets [3, H^, [2l| , before crossing over to a 
simple aging form, Tq, ~ tw, at large times, tv/ ^ Ta- In 
the initial regime, we also found that the energy density 
decays very quickly, suggesting that the system evolves 
rapidly from its highly disordered initial state to form 
a nearly frozen structure, which then ages more slowly. 
Thus, this initial regime is likely quite sensitive to the 
detailed initialization procedure of the system both in 
simulations and experiments. The 'universal' features of 
the aging dynamics of concentrated hard spheres should 
be discussed for the second regime only, and our simu- 
lations indicate that concentrated hard spheres follow a 
simple aging form, as indeed found for very many glassy 
materials [l|. 

In the aging regime, we found that particle motions 
are very restrained, with particles moving on average less 
than their own diameter over the entire duration of sim- 
ulation. In the aging regime, ^ r^, wc found that 
the self-intermediate functions decay algebraically at long 
times, in agreement with several simulations of soft and 
hard particle systems [1^, [13 • This is also in agreement 
with recent experiments performed on colloidal suspen- 
sions that are dense enough, so that the aging regime 
does not cross over at long times towards equilibrium 
behavior 

We also found that single particle motion is sub- 
diffusive, Ar'^{tyf,T) ~ r'^, with v < 1. This behavior 
is actually quite consistent with the results obtained by 
optical and confocal microsopy on hard and soft colloids, 
where diffusive behavior is not reached in the experimen- 
tal time window [l3, [H, [13] ■ This result is also in broad 
agreement with earlier numerical work [25| , although sub- 
diffusion was never described in detail before. Wc have 
suggested that a natural theoretical framework to ana- 
lyze our results should the Aging Continuous Time Ran- 
dom Walk (ACTRW), i.e. the nonequilibrium extension 
of the random walk picture used to describe single parti- 
cle motion at thermal equilibrium near the glass transi- 
tion 0, [31 . This formalism makes a number of detailed 
predictions for the aging regime. Here, for lack of statis- 
tics, we have simply been able to show that our data 
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are in qualitatively good agreement with this picture, 
which can link the shape of the distribution of parti- 
cle displacements to the subdifFusion exponent. It would 
be extremely interesting to extend this analysis in future 
work and check in more detail how far the ACTRW pic- 
ture can be pushed to describe the aging dynamics of 
concentrated hard spheres. 

We finally showed, by resolving the measurement of 
time correlation functions in space and time, that the 
aging dynamics of concentrated hard spheres is highly in- 
termittent, with very sudden relaxation events separating 
long periods of time where very little dynamics occurs. 
Similar dynamic events have been observed in numeri- 
cal work before, and were coined 'earthquakes' [IHl or 
'avalanches' [s^. Although these names and additional 
numerical evidence s ugg ested that these events could well 
be system-spanning [33| , we showed by using much larger 
system sizes that these events remained actually spatially 
localized and do not grow with system size. This is in 
fact quite consistent with the experimental report of dy- 
namic clusters that grow rather modestly in aging col- 
loidal samples [isl [T7| , and with numerical work report- 
ing a slow growth of four-point dynamic susceptibilities 
in aging Lennard- Jones glasses [3l| . 

Thus, we find that in concentrated hard spheres the dy- 



namics are intermittent as it occurs in several more com- 
plex colloidal materials, but this leads neither to anoma- 
lous compressed exp onential relaxation for time corre- 
lation functions Il4ll . nor to ballistic motion over large 
distances [13, [l6l. llSl [Tgl . [53 | , whose origin thus remains 
largely mysterious. This means that model systems such 
as hard spheres or Lennard-Jones glasses are not good 
starting points to gain insight into the nature of these in- 
triguing aging dynamics. It is not clear what ingredient 
these models are missing, since similar anomalous aging 
dynamics was recently reported in simple 'hard' molecu- 
lar glasses approaching the glass transition [2^, [2^. It 
would be very interesting to discover a model system 
displaying the same type of anomalous aging dynamics, 
with a particle-scale dynamics that can be followed in the 
simulations in the way similar to what we did for hard 
spheres. 
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